R Notebook

#SF Mapping of Edmonton. 
#Install packages, library

# install.packages("sf")
# install.packages("tmap")
# install.packages("ggmap")
# install.packages("ggthemes")
# install.packages("rgeos")
# install.packages("grid")
# install.packages("gridExtra")
# install.packages("reshape2")
# install.packages("scales")
# install.packages(c("rgdal", "htmlwidgets"), dependencies = TRUE)
#install.packages("rgdal")
#install.packages("geojsonsf")


library(geojsonsf)

package 㤼㸱geojsonsf㤼㸲 was built under R version 4.0.5

library(sf)
library(sp)
library(tmap)

package 㤼㸱tmap㤼㸲 was built under R version 4.0.5

library(dplyr)
library(maptools)

Checking rgeos availability: TRUE

library(ggmap)

package 㤼㸱ggmap㤼㸲 was built under R version 4.0.5Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
Please cite ggmap if you use it! See citation("ggmap") for details.

library(forcats)
library(ggthemes)

package 㤼㸱ggthemes㤼㸲 was built under R version 4.0.5

library(rgeos)

package 㤼㸱rgeos㤼㸲 was built under R version 4.0.5rgeos version: 0.5-5, (SVN revision 640)
 GEOS runtime version: 3.8.0-CAPI-1.13.1 
 Linking to sp version: 1.4-5 
 Polygon checking: TRUE 

library(broom)
library(plyr)

--------------------------------------------------------------------------------------------------------------------------------------------------
You have loaded plyr after dplyr - this is likely to cause problems.
If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
library(plyr); library(dplyr)
--------------------------------------------------------------------------------------------------------------------------------------------------

Attaching package: 㤼㸱plyr㤼㸲

The following objects are masked from 㤼㸱package:dplyr㤼㸲:

    arrange, count, desc, failwith, id, mutate, rename, summarise, summarize

library(grid)
library(gridExtra)

package 㤼㸱gridExtra㤼㸲 was built under R version 4.0.5
Attaching package: 㤼㸱gridExtra㤼㸲

The following object is masked from 㤼㸱package:dplyr㤼㸲:

    combine

library(reshape2)

package 㤼㸱reshape2㤼㸲 was built under R version 4.0.5
Attaching package: 㤼㸱reshape2㤼㸲

The following objects are masked from 㤼㸱package:data.table㤼㸲:

    dcast, melt

library(scales)
library("rgdal")

package 㤼㸱rgdal㤼㸲 was built under R version 4.0.5rgdal: version: 1.5-23, (SVN revision 1121)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
Path to GDAL shared files: G:/Documents/R/win-library/4.0/rgdal/gdal
GDAL binary built with GEOS: TRUE 
Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
Path to PROJ shared files: G:/Documents/R/win-library/4.0/rgdal/proj
PROJ CDN enabled: FALSE
Linking to sp version:1.4-5
To mute warnings of possible GDAL/OSR exportToProj4() degradation,
use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
Overwritten PROJ_LIB was G:/Documents/R/win-library/4.0/rgdal/proj

library(tmap)
library(leaflet) 

package 㤼㸱leaflet㤼㸲 was built under R version 4.0.5

#use net map
 
#manually definen bounding box of edmonton
edm_map <- get_map(location =  c(left = -113.702 , bottom = 53.366 , right = -113.317, top = 53.6855 ), 
        zoom = 10,
        scale = 12,
        maptype = "toner",
        source = "osm"
        )

ggmap(edm_map)

#set up map parameters
mapTheme <- function(base_size = 12) {
  theme(
    text = element_text( color = "black"),
    plot.title = element_text(size = 18,colour = "black"),
    plot.subtitle=element_text(face="italic"),
    plot.caption=element_text(hjust=0),
    axis.ticks = element_blank(),
    panel.background = element_blank(),
    panel.grid.major = element_line("grey80", size = 0.1),
    strip.text = element_text(size=12),
    axis.title = element_blank(),
    axis.text = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.grid.minor = element_blank(),
    strip.background = element_rect(fill = "grey80", color = "white"),
    plot.background = element_blank(),
    legend.background = element_blank(),
    legend.title = element_text(colour = "black", face = "italic"),
    legend.text = element_text(colour = "black", face = "italic"))
}
# Define some palettes
palette_9_colors <- c("#0DA3A0","#2999A9","#458FB2","#6285BB","#7E7CC4","#9A72CD","#B768D6","#D35EDF","#F055E9")
palette_8_colors <- c("#0DA3A0","#2D97AA","#4D8CB4","#6E81BF","#8E76C9","#AF6BD4","#CF60DE","#F055E9")
palette_7_colors <- c("#2D97AA","#4D8CB4","#6E81BF","#8E76C9","#AF6BD4","#CF60DE","#F055E9")
palette_1_colors <- c("#0DA3A0")
 
# Map it
bmMap <- ggmap(edm_map) + mapTheme() + 
  labs(title="Edmonton basemap")
bmMap

NA
NA
NA
NA

#Create map of accessibility based on data crunched by R5R Edmonton TTM combined with neighborhood shapefile.

#Explanation: https://www.infoworld.com/article/3505897/how-to-do-spatial-analysis-in-r-with-sf.html
#sf documentation: https://cran.r-project.org/web/packages/sf/sf.pdf

#Read in as poly file
Edm_hoods <- st_read("G:/Desktop/Folders/Edu/Summer 2021/Data Files/ArcGIS Files/School Data/COUNT_ELEMENTARY_schools_per_hood.shp")

Reading layer `COUNT_ELEMENTARY_schools_per_hood' from data source `G:\Desktop\Folders\Edu\Summer 2021\Data Files\ArcGIS Files\School Data\COUNT_ELEMENTARY_schools_per_hood.shp' using driver `ESRI Shapefile'
Simple feature collection with 400 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 18949.21 ymin: 5911407 xmax: 48145.16 ymax: 5953685
Projected CRS: NAD27 / Alberta 3TM ref merid 114 W

plot(Edm_hoods) #plots all columns on the map. Can show number of schools per neighborhood

# Define the bounding box
# bbox <- Edm_hoods@bbox

neighborhood_access <- acc_20min #from r5 notebook
names(neighborhood_access) <- c("id", "accessible")

Edm_hoods_dt <- as.data.table(Edm_hoods)
names(Edm_hoods_dt) <- c("id","name","number","lon","lat","count","geometry")
Edm_hoods_dt$id <- as.character(Edm_hoods_dt$id)#need to change id to char.


Edm_hoods_dt <- left_join(Edm_hoods_dt, neighborhood_access, by = "id") #aggregated and joined data set
Edm_hoods_dt <- as.data.table(Edm_hoods_dt)
Edm_hoods_dt

Edm_hoods_dt_sf <- st_as_sf(Edm_hoods_dt) #turns it back to SF object
plot(Edm_hoods_dt_sf["accessible"], key.pos = 1) #plot the count column. Number of schools accessible within 20 minutes

NA
NA
NA

#interactive map

#resource https://cengel.github.io/R-spatial/mapping.html
#tmap documentation: https://rdrr.io/cran/tmap/man/tm_polygons.html

tm_shape(Edm_hoods_dt_sf) +
  tm_polygons("accessible", 
              style="pretty", #or "cat" for actual count 
              title="Number of public\nschools accessible\nin 20 mins by neighborhood",
              interactive=TRUE,
              popup.vars="name",
              popup.format=list()
  )

tmap_mode("plot") #changes between view mode (interactive) and plot mode (static)

tmap mode set to plotting

#Interactive mapping with leaflet
#can be used to create mobile friendly java based maps.
#reference: https://cengel.github.io/R-spatial/mapping.html

#need to manually repreject the map to get it to match with web projection
Edm_hoods_WGS84 <- st_transform(Edm_hoods_dt_sf, 4326)

leaflet(Edm_hoods_WGS84) %>%
  addPolygons()

pal_fun <- colorBin("YlOrRd", NULL, bins = 5) #colour options

p_popup <- paste0(Edm_hoods_WGS84$name, "<strong> - Schools reachable in 20 mins: </strong>", Edm_hoods_WGS84$accessible) #popup options

interactive_edmSchool_map <- leaflet(Edm_hoods_WGS84) %>%
  addPolygons(
    stroke = FALSE, # remove polygon borders
    fillColor = ~pal_fun(accessible), # set fill color with function from above and value
    fillOpacity = 0.6, smoothFactor = 0.5, # make it nicer
    popup = p_popup, 
    group = "Edm school data") %>% # add popup
  
  addTiles() %>% #automatically add basemap from OSM
  addLegend("bottomright",  # location
            pal=pal_fun,    # palette function
            values=~accessible,  # value to be passed to palette function
            title = 'Number of public\nschools accessible\nin 20 mins by neighborhood') %>% # legend title

  addLayersControl(baseGroups = c("OSM", "Carto"), #turn on or off data layer
                   overlayGroups = c("Edm school data"))  


interactive_edmSchool_map

NA

st_write(Edm_hoods_dt, "G:/Desktop/Folders/Edu/Summer 2021/Data Files/R Code/Edm_HighSchool_Access.shp") #writes the sf file to csv

Writing layer `Edm_HighSchool_Access' to data source `G:/Desktop/Folders/Edu/Summer 2021/Data Files/R Code/Edm_HighSchool_Access.shp' using driver `ESRI Shapefile'
Writing 400 features with 7 fields and geometry type Multi Polygon.